Stiffness characteristics analysis of a Biglide industrial parallel robot considering the gravity of mobile platform and links

For the machining process of industrial parallel robots, the gravity generated by the weight of mobile platform and links will lead to the deviation of the expected machining trajectory of the tool head. In order to evaluate this deviation and then circumvent it, it is necessary to perform the robotic stiffness model. However, the influence of gravity is seldom considered in the previous stiffness analysis. This paper presents an effective stiffness modeling method for industrial parallel robots considering the link/joint compliance, the mobile platform/link gravity, and the mass center position of each link. First, the external gravity corresponding to each component is determined by the static model under the influence of gravity and mass center position. Then, the corresponding Jacobian matrix of each component is obtained by the kinematic model. Subsequently, the compliance of each component is obtained by cantilever beam theory and FEA-based virtual experiments. In turn, the stiffness model of the whole parallel robot is determined and the Cartesian stiffness matrix of the parallel robot is calculated at several positions. Moreover, the principal stiffness distribution of the tool head in each direction over the main workspace is predicted. Finally, the validity of the stiffness model with gravity is experimentally proved by the comparison of the calculated stiffness and measured stiffness in identical conditions.

For the machining process of industrial parallel robots, the gravity generated by the weight of mobile platform and links will lead to the deviation of the expected machining trajectory of the tool head. In order to evaluate this deviation and then circumvent it, it is necessary to perform the robotic stiffness model. However, the influence of gravity is seldom considered in the previous stiffness analysis. This paper presents an effective stiffness modeling method for industrial parallel robots considering the link/joint compliance, the mobile platform/link gravity, and the mass center position of each link. First, the external gravity corresponding to each component is determined by the static model under the influence of gravity and mass center position. Then, the corresponding Jacobian matrix of each component is obtained by the kinematic model. Subsequently, the compliance of each component is obtained by cantilever beam theory and FEA-based virtual experiments. In turn, the stiffness model of the whole parallel robot is determined and the Cartesian stiffness matrix of the parallel robot is calculated at several positions. Moreover, the principal stiffness distribution of the tool head in each direction over the main workspace is predicted. Finally, the validity of the stiffness model with gravity is experimentally proved by the comparison of the calculated stiffness and measured stiffness in identical conditions. In many modern applications, such as manufacturing robots 1,2 , bionic robots 3,4 , medical robots 5,6 and aeronautical robots 7,8 , the mechanism is subject to external payloads and gravity of its own components which in turn induce deformations causing positioning errors. In order to improve motion accuracy and machining accuracy, stiffness analysis becomes an important solution. Moreover, stiffness analysis is also very important for the design stage of a parallel robot 9 . This is why stiffness analysis has become a research hotspot in recent years.
In terms of stiffness modeling, stiffness modeling methods can be mainly divided into three categories 10 : the finite element analysis method (FEA), the matrix structural analysis method (MSA) and the virtual joint modeling method (VJM). The FEA divides the target model into smaller and simpler connected element modules according to its real dimension and shape. Therefore, the modeling accuracy of FEA is very high, but its calculation cost is also very large. It needs the help of commercial FEA software to complete the modeling of the target model, and in order to ensure that the nodes of the module in the division match the nodes on the model, it needs to re-mesh the finite element model repeatedly. That's why it's expensive to calculate, and the process is time-consuming and tedious. Therefore, the FEA is not suitable for parametric stiffness model, which requires stiffness evaluation of the entire workspace.
The MSA method takes each component of the robot as a structural unit, which is a large compliant element, and combines the main ideas of the FEA to establish stiffness model of the robot. The main idea of this method is to obtain the stiffness matrix of the robot by combining the stiffness matrix of standard elements such as links, joints and frames with matrix analysis method. However, the distortion of this method is that these substructures are regarded as regular shapes, but the stiffness matrix of their standard element cannot accurately describe the actual stiffness of these substructures. Although some scholars have made some innovations on this method [11][12][13] , this method is still not suitable for the stiffness modeling of robot with the substructure of multi-node element.
In these approaches, the VJM is the most attractive and is adopted in this paper. This method obtains the stiffness model of the robot by establishing the force and deformation mapping relationship between joint space and operation space and applying Jacobian matrix and Hooke's law. And the overall compliance factor comes www.nature.com/scientificreports/ 3-RR parallelogram-link structure (j = 1, 2, 3). For accuracy of expression, i and j are taken as such throughout the whole paper. The three links in the parallelogram-link structure are connected in parallel, then the Cartesian stiffness matrix K i C,Pl of the 3-RR parallelogram-link structure can be expressed as where K i θ,Pl is the stiffness matrix of the 3-RR parallelogram-link structure, K i G,Pl denotes the stiffness matrix due to the gravity of the mobile platform for the 3-RR parallelogram-link structure, J i θ,Pl is the Jacobian matrix of the virtual spring corresponding to the 3-RR parallelogram-link structure.
(2) Stiffness model of the single limb (P(3-RR) mechanism) In the P(3-RR) mechanism, prismatic joint (actuator), slidable platform (SP) and parallelogram-link structure are connected in series, then the Cartesian stiffness matrix K i C,limb of the single limb can be expressed as where K i C,Ac is the Cartesian stiffness matrix of the actuator (lead screw driving system), K i C,SP is the Cartesian stiffness matrix of the slidable platform, K θ,Ac is the stiffness coefficient of the driving system consisting of the motor and the screw, K i θ,SP is the stiffness matrix of the slidable platform, K i G,SP denotes the stiffness matrix due to the gravity of the mobile platform and links for the slidable platform, J i θ,Ac is the Jacobian matrix related to the actuator (prismatic joint), J i θ,SP is the Jacobian matrix of the virtual spring corresponding to the slidable platform.
(3) Stiffness model of the parallel limbs (2-P(3-RR) mechanism) In the 2-P(3-RR) mechanism, the two limbs are connected in parallel and the size of the mobile platform is considered, then the Cartesian stiffness matrix K C,limb of the 2P(3-RR) mechanism can be expressed as where J i v defines geometrical mapping between end-points of serial limbs and reference point frame (end-effector). (4) Stiffness model of the whole Biglide parallel robot (P[2-P(3-RR)] mechanism). www.nature.com/scientificreports/ In the P[2-P(3-RR)] mechanism, prismatic joint (actuator) and 2-P(3-RR) mechanism are connected in series in the Biglide parallel robot, then the Cartesian stiffness matrix K C of the Biglide parallel robot can be expressed as where J 0 θ,Ac is the Jacobian matrix corresponding to the first prismatic joint. Moreover, in order to solve the above equations, the static and kinematic analysis of the Biglide parallel robot is required.

Static analysis of the Biglide parallel robot
In this section, the respective static models are developed at each component of the Biglide parallel robot, taking into account the gravity generated by the weight of the mobile platform and the moveable links [26][27][28][29] . The mathematical equation of the static force/torque corresponding to each component is derived, and thus the gravity received at the component is obtained.
Coordinate system of the Biglide parallel robot. The mechanical structure illustration of the Biglide parallel robot is offered in Fig. 1. The Biglide parallel robot is a two degree-of-freedom (2-dof) parallel mechanism (PM) and consist of two modular parallelogram-links connected to the mobile platform (MP) and two sliding units which are actuated via lead screw system, where the sliding units are installed on the upper fixed platform. The Biglide parallel robot can generate two translations in the horizontal/vertical directions within a single plane, and the single moving plane can also move along its normal direction.
To convenience the calculation, a sequence of coordinate systems is attached to each substructure of the Biglide parallel robot. where R(.) represents a 3 × 3 rotation matrix that rotates the corresponding angle around the corresponding axis, R O is an identity matrix, and R n is the orientation matrix of the coordinate system {O n } with respect to the coordinate system {O}.

Static equilibrium equation of points M i .
In order to evaluate the influence of the mobile platform weight at points M 1 and M 2 , it is necessary to analysis the static force/torque balance at points M 1 and M 2 , as shown in Fig. 3. The static equilibrium equations at point M can be displayed as  www.nature.com/scientificreports/ ρ Bi1 , ρ Bi2 and ρ Bi3 denote the reaction wrench applied at joints B i1 , B i2 and B i3 , respectively; η Bi1 , η Bi2 and η Bi3 denote the equivalent external wrench for joints B i1 , B i2 and B i3 applied at point M, respectively; s(r Bi1 )/s(r Bi2 ) /s(r Bi3 ) represents the skew-matrix of the position vector r Bi1 /r Bi2 /r Bi3 .

Static equilibrium equation of joints A 1j .
In order to evaluate the influence of the weight of the mobile platform and links on the joint A ij , it is necessary to analysis the static force/torque balance at point A ij , as shown in Fig. 3b. The static equilibrium equations at point B ij can be displayed as where m i1 , m i2 and m i3 are the mass of links L i1 , L i2 and L i3 , respectively; G i1 , G i2 and G i3 are the mass center of links L i1 , L i2 and L i3 , respectively; f Ai1 and τ Ai1 denote the reaction force and torque applied at point A i1 , f Ai2 and τ Ai2 denote the reaction force and torque applied at point A i2 , f Ai3 and τ Ai3 denote the reaction force and torque applied at point A i3 ; r Bi1 Gi1 denotes the position vector from point B i1 to point G i1 expressed in the coordinate system {B i1 }, r Bi2 Gi2 denotes the position vector from point B i2 to point G i2 expressed in the coordinate system {B i2 }, r Bi3 Gi3 denotes the position vector from point B i3 to point G i3 expressed in the coordinate system {B i3 }; r Bi1 Ai1 denotes the position vector from point B i1 to point A i1 expressed in the coordinate system {B i1 }, r Bi2 Ai2 denotes the position vector from point B i2 to point A i2 expressed in the coordinate system {B i2 }, r Bi3 Ai3 denotes the position vector from point B i3 to point A i3 expressed in the coordinate system {B i3 }.
According to the decoupling of Eq. (11), the reaction force acting on the joints A i1 , A i2 and A i3 can be obtained where η Ai1 , η Ai2 and η Ai3 denote the resultant external wrench for joint A i1 , A i2 and A i3 applied at point M, respectively; η i1 , η i2 and η i3 denote the equivalent gravitational wrench of links L i1 , L i2 and L i3 applied at point M, www.nature.com/scientificreports/ respectively; ρ Ai1 , ρ Ai2 and ρ Ai3 denote the reaction wrench applied at joints A i1 , A i2 and A i3 , respectively; s(r Ai1 ) /s(r Ai2 )/s(r Ai3 ) represents the skew-matrix of the position vector r Ai1 /r Ai2 /r Ai3 .

Kinematic analysis of the Biglide parallel robot
In this section, the kinematic analysis of the Biglide parallel robot is carried out so that the Jacobian matrix of each substructure is derived [27][28][29][30][31][32][33] . Due to the position relationship and series/parallel relationship of each substructure, the kinematic equilibrium equation is established for each substructure of the Biglide parallel robot.

Kinematic modeling for 3-RR parallelogram-link structure.
In order to obtain the Jacobian matrices of each kinematic chain in the 3-RR, it is necessary to analyze the kinematic transitions from point A i to B i through each chain, as shown in Fig. 3b. The kinematic homogenous matrix of the first chain in the left 3-RR Pl can be expressed as where P(.) represents a 3 × 1 displacement vector that translates the corresponding distance along the corresponding axis, I 3 is a 3 × 3 identity matrix.
The position and direction of the end point B1 can be extracted from the matrix T 1 Pl,1 in a standard way 27,33,34 , whereby the kinematic model can be rewritten in the form of a vector function as follows  /s ω 1 q,Pl,1k is represents the skew-matrix of the orientation vector ω 1 θ,Pl,1n /ω 1 q,Pl,1k . According to the above theory, the Jacobian matrices of the first chain in the left 3-RR Pl can be obtained by deriving Eq. (15) as follows: where S q and C q represent sine and cosine of the angle q 1 , respectively.
Similarly, the Jacobian matrices for the second and third chains in the left 3-RR Pl, respectively, can be expressed as For the right 3-RR Pl, the Jacobian matrices of its three chains can also be obtained in the same way.
Kinematic modeling for the single limb (P(3-RR) mechanism). In order to obtain the Jacobian matrices of each kinematic limb, it is necessary to analyze the kinematic transitions from point O i to M i through each limb, as shown in Fig. 3a ; where p 1 θ,limb,n denote the position vector associated with θ 1 n ′ , s ω 1 θ,limb,n is represents the skew-matrix of the orientation vector ω 1 θ,limb,n . According to the above theory, the Jacobian matrices of limb 1 can be obtained by deriving Eq. (21) as follows: For limb 2, the Jacobian matrices of its three components can also be obtained in the same way.

Kinematic modeling for the whole robot (P[2-P(3-RR)] mechanism). The Biglide parallel robot is
composed of a lead screw system (prismatic joint) capable of producing translation along the y O -axis in series with a 2-P(3-RR) parallel structure, where the 2-P(3-RR) parallel structure includes two limbs, a screw system capable of producing translation along the x O -axis, and a mobile platform, as shown in Fig. 2. In order to calculate the stiffness of the whole robot, the kinematic Jacobian matrices of the mobile platform and the prismatic joint need to be evaluated. These Jacobian matrices can be expressed as

Numerical analysis
In this section, the principal stiffness of the Biglide parallel robot considering the gravity of the mobile platform and the links is investigated numerically using the developed stiffness model. The dimensional and mass parameters of the robot are shown in Table 1. The position parameters of each limb are shown in Tables 2 and 3. The dimensional parameters of the link are shown in Table 4. The reachable workspace of Biglide parallel robot in the x-z plane is the yellow area in Fig. 4, and the blue area enclosed by W 1 -W 4 is its main workspace, where Q 1 and Q 2 are the midpoints of W 1 W 2 and W 3 W 4 , respectively.
Compliance matrices. The 3-RR Pl contains three identical cylindrical links, as shown in Fig. 3. In order to calculate the Cartesian stiffness matrix of the 3-RR Pl, it is necessary to obtain the compliance matrix of each link by means of the cantilever beam principle 16,19 . Its compliance matrix K i −1 θ,Pl,j can be expressed as where I y , I z and I p are the quadratic and polar moments of inertia of the cross-section, and E and G are the Young's and shear modules, respectively. Moreover, the slidable platform stiffness can be identified on the basis of CAD models and FEA methods 10,25,37 . Its compliance matrix K i −1 θ,SP can be expressed as (unit: N, m, rad)   Distribution of the principal stiffness. In order to better evaluate the Cartesian stiffness of the Biglide parallel robot, we visualize the distribution of its principal stiffness in six directions in its workspace, as shown in Figs. 5 and 6. K tx , K ty , K tz and K rx , K ry , K rz are the diagonal elements of the Cartesian stiffness matrix K O C , which represent the principal translational stiffness along the x, y, z axes and the principal rotational stiffness around the x, y, z axes, respectively. Figure 5a shows the principal translational stiffness K tx along the x-axis with and without gravity, it can be seen that the stiffness with gravity improves compared to the stiffness without gravity, and they both decrease as the z-axis value increases; Fig. 5b shows the principal translational stiffness K ty along the y-axis with and without gravity, it is obvious that the stiffness with gravity is the same as the stiffness without gravity, and this stiffness is a constant; Fig. 5c shows the principal translational stiffness K tz along the z-axis with and without gravity, it has similar stiffness characteristics to K tx ; and Fig. 5d shows the mean difference between the principal translational stiffness with and without gravity. Moreover, among the principal translational stiffness K tz is the highest order of magnitude translational stiffness, K tx is the second and K ty is the smallest. Figure 6 is similar to Fig. 5, but the difference is that Fig. 6 shows the principal rotational stiffness with and without gravity and their mean difference. It is worth noting that the stiffness K rx in Fig. 6a becomes smaller due to the gravitational influence and it is the only principal stiffness where gravity has a negative effect. K rx increases with increasing z-axis values, while K ry and K rz decrease with increasing z-axis values. Furthermore,   www.nature.com/scientificreports/ K ry is the highest order of magnitude rotational stiffness in the principal rotational stiffness, and K rx and K rz have the same order of magnitude.

Experimental analysis
To validate the correctness of the stiffness model with gravity, the knocking experiments were employed on the Biglide industrial parallel robot. Using a laser Doppler vibrometer as the main experimental equipment, the principal stiffness of the mobile platform center of the Biglide parallel robot was measured in all directions, as shown in Fig. 7. Table 5 lists the principal stiffness in each direction at point Q 1 . The values at point Q 2 are listed in Table 6. These values are presented in a more intuitive form in Fig. 8. It can be found that the principal stiffness of the stiffness model with gravity is very close to the experimental results. The root mean square (RMS) values between the developed model and the experimental results at points Q 1 and Q 2 are 2.67 and 2.77%, respectively 38,39 . And the RMS values between the stiffness model without gravity and the experimental results at Q 1 and Q 2 are 9.043 and 9.353%, respectively. The experimental results show a good agreement with the theoretical values of the 'with gravity' model. Therefore, it can be concluded that the stiffness model with gravity is effective, which lays a solid foundation for the real applications.

Conclusions
This paper proposes a novel stiffness modeling method for an industrial parallel robot, which is important for the machining of parallel robots. The main research results are summarized as follows: (1) The influences of the gravity of the movable links, the gravity of the mobile platform and their corresponding mass centers are considered simultaneously. In order to obtain the values of gravity, the mathematical equations of the static model were established for each joint. (2) The stiffness model of each component is established by considering its compliance matrix, the compliance of the corresponding joint, its Jacobian matrix from the corresponding mathematical equations of the kinematic model, and the corresponding stiffness matrix due to external gravity. For the accuracy of the calculation, the compliance matrix of the irregular component is identified by using the FEA-based virtual experiment.  www.nature.com/scientificreports/  www.nature.com/scientificreports/ (3) The validity of the stiffness model with gravity was verified by comparing the calculated stiffness of the Biglide industrial parallel robot with the knockout experiment. The results show that the method is able to obtain sufficient computational accuracy to predict the stiffness distribution in the task workspace. This enables the robot to obtain more accurate tool head trajectory before machining, thus improving its machining accuracy. (4) The stiffness modeling method is also applicable to the stiffness estimation of the robot structure design stage, which provides effective help for the improvement and verification of the structure.

Data availability
The datasets generated and/or analysed during the current study are not publicly available due the data also forms part of an ongoing study but are available from the corresponding author on reasonable request.